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Abstract.  A  weak  form  Galerkin  finite  element  model  for  the  nonlinear  quasi-static  and 
fully  transient  analysis  of  initially  straight  viscoelastic  beams  is  developed  using  the  kine¬ 
matic  assumptions  of  the  third-order  Reddy  beam  theory.  The  formulation  assumes  linear 
viscoelastic  material  properties  and  is  applicable  to  problems  involving  small  strains  and 
moderate  rotations.  The  viscoelastic  constitutive  equations  are  efficiently  discretized  using 
the  trapezoidal  rule  in  conjunction  with  a  two-point  recurrence  formula.  Locking  is  avoided 
through  the  use  of  standard  low  order  reduced  integration  elements  as  well  through  the 
employment  of  a  family  of  elements  constructed  using  high  polynomial-order  Lagrange  and 
Hermite  interpolation  functions. 


1.  Introduction 

Materials  exhibiting  characteristics  of  both  elastic  solids  as  well  as  viscous  fluids  are  com¬ 
monly  known  as  viscoelastic  materials.  Prominent  examples  include  metals  at  elevated  tem¬ 
peratures,  polymers,  rubbers  and  concrete.  The  theoretical  foundations  of  viscoelasticity 
are  well  established.  We  refer  to  the  standard  texts  of  Fliigge  [14],  Christensen  [9],  Findley 
[13]  and  Reddy  [36]  for  an  overview  on  the  theory  of  viscoelastic  material  behavior,  as  well 
as  the  classical  analytical  solution  techniques  that  may  be  used  to  solve  simple  viscoelastic 
boundary  value  problems.  Viscoelastic  materials  are  often  highly  desirable  for  use  in  struc¬ 
tural  components,  clue  to  their  natural  ability  to  dampen  out  structural  vibrations.  The 
capability  to  satisfactory  predict  the  mechanical  response  of  viscoelastic  structures  therefore 
becomes  of  great  importance  in  engineering  design  scenarios. 

A  variety  of  beam  theory  based  finite  element  models  have  been  presented  in  the  literature 
for  the  analysis  of  viscoelastic  structures.  The  majority  of  these  formulations  employ  some 
form  of  either  the  Euler-Bernoulli  or  Timoshenko  beam  theories  and  are  mostly  restricted 
to  small  strain  analysis.  The  formulations  differ  in  how  the  convolution  form  of  the  vis¬ 
coelastic  constitutive  equations  are  temporally  discretized.  A  popular  approach  adopted  by 
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many  researchers  is  to  employ  the  Laplace  transform  method  directly  in  the  construction 
of  the  finite  element  equations  [8,  1,  39].  In  this  approach,  terms  associated  with  the  time 
domain,  including  the  convolution  integral,  are  transformed  into  variables  associated  with 
the  Laplace  space  s.  A  successful  numerical  simulation  therefore  requires  an  efficient  and 
accurate  inversion  of  the  solution  in  s  space  back  to  the  time  domain.  Many  of  the  key  ideas 
are  presented  in  work  of  Akoz  and  Kadioglu  [1],  wherein  a  Timoshenko  beam  element  is  de¬ 
veloped  using  mixed  variational  principles.  In  their  work,  the  finite  element  model  requires 
numerical  inversion  from  the  Laplace-Carson  domain  back  to  the  time  domain.  Temel  et  al. 
[39]  utilized  the  Durbin’s  inverse  Laplace  transform  method  in  their  analysis  of  cylindrical 
helical  rods  (based  on  the  Timoshenko  beam  hypotheses). 

Additional  numerical  formulations  for  viscoelastic  beams  are  based  on  the  Fourier  trans¬ 
form  method  [7],  the  anelastic  displacement  (ADN)  procedure  [40,  28]  and  the  Golla-Hughes- 
McTavish  (GHM)  method  [24,  25,  4,  3].  It  has  been  noted  that  when  the  relaxation  moduli 
are  given  as  Prony  series,  the  convolution  form  of  linear  viscoelastic  constitutive  equations 
may  be  equivalently  expressed  as  a  set  of  ordinary  differential  equations  in  terms  of  a  col¬ 
lection  of  internal  strain  variables.  Numerical  discretization  procedures  exploiting  this  ODE 
form  of  the  viscoelastic  constitutive  equations  have  been  successfully  adopted  in  the  works 
of  Johnson  et  al.  [19]  and  Austin  and  Inman  [2], 

The  formulations  described  above  are  restricted  to  a  class  of  problems  involving  infin¬ 
itesimal  strains  and  small  deflections.  Among  the  finite  element  formulations  appearing 
in  the  literature  for  nonlinear  viscoelastic  shell  structures  are  the  works  of  Kennedy  [21], 
Oliveira  and  Creus  [27]  and  Hammerand  and  Kapania  [16].  More  recently,  Payette  and 
Reddy  [29]  presented  quasi-static  finite  element  formulations  for  Euler-Bernoulli  and  Timo¬ 
shenko  beams  that  allow  for  loading  scenarios  that  produce  large  transverse  displacements, 
moderate  rotations  and  small  strains.  This  previous  work  may  be  viewed  as  a  bridge  between 
the  formulations  associated  with  either:  (a)  small  strains  or  (b)  fully  finite  deformations.  In 
this  prior  work,  the  trapezoidal  rule  was  employed  in  conjunction  with  a  two-point  recur¬ 
rence  formula  for  efficient  temporal  integration  of  the  viscoelastic  constitutive  equations. 
The  objective  of  the  present  paper  is  to  extend  the  work  of  Payette  and  Reddy  to  create 
an  efficient  locking-free  nonlinear  finite  element  framework  for  the  analysis  of  viscoelastic 
beam  structures  based  on  the  high-order  Reddy  beam  theory  [30,  18,  41]  for  use  in  both 
quasi-static  as  well  as  fully  dynamic  analysis. 
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2.  Kinematics  of  deformation 


2.1.  The  displacement  field.  There  are  a  variety  of  beam  theories  that  have  been  success¬ 
fully  employed  in  the  mechanical  analysis  of  structural  elements.  Such  theories  are  typically 
formulated  in  terms  of  truncated  Taylor  series  expansions  of  the  components  of  the  dis¬ 
placement  field  with  respect  to  the  thickness  coordinate.  The  most  simple  and  commonly 
employed  theories  are  the  Eulcr-Bernoulli  beam  theory  (EBT)  and  the  Timoshenko  beam 
theory  (TBT).  The  major  deficiency  associated  with  the  EBT  is  failure  to  account  for  defor¬ 
mations  associated  with  shearing.  The  TBT  relaxes  the  normality  assumption  of  the  EBT 
and  admits  a  constant  state  of  shear  strain  on  a  given  cross-section.  As  a  result,  the  TBT 
necessitates  the  use  of  shear  correction  coefficients  in  order  to  accurately  predict  transverse 
displacements.  The  third-order  Reddy  beam  theory  (RBT)  was  introduced  to  both  account 
for  the  effects  of  shear  strains  and  to  also  produce  a  parabolic  variation  of  the  shear  strain 
through  the  thickness  [30,  18,  41].  As  a  result,  in  the  RBT  there  is  no  need  to  introduce 
shear  correction  coefficients. 

Before  presenting  the  displacement  field  associated  with  the  RBT  we  first  introduce  some 
standard  notation.  We  let  B  C  M3,  an  open  and  bounded  set,  denote  the  material  or 
reference  configuration  of  the  beam.  The  material  configuration  may  be  expressed  as  B  = 
SlxA,  where  hi  =  (0,  L)  and  L  is  the  initial  length  of  the  beam.  In  addition  the  quantity 
A  represents  the  beam’s  cross-sectional  area.  A  typical  material  point  belonging  to  B  is 
denoted  as  X  =  (X,Y,Z).  Likewise  the  spatial  or  current  configuration  of  the  beam  is 
denoted  by  Bt  and  the  associated  points  are  expressed  as  x  =  (x,  y,  z).  Points  in  the  spatial 
configuration  are  related  to  points  in  the  material  configuration  by  the  standard  bijective 
mapping  %  :  B  x  M  — >  Bt.  As  a  result  x  =  %(X,  t).  The  displacement  field  associated  with 
the  mapping  may  be  expressed  in  the  usual  manner  as  u(X,f)  =  %(X,  t)  —  X. 

In  the  present  work,  we  constrain  the  displacement  field  to  conform  to  the  kinematic 
assumptions  of  the  Reddy  beam  theory.  The  displacement  field  in  the  Reddy  beam  theory 
(for  a  beam  with  a  rectangular  cross  section)  is  taken  as 

u(X,  Y,  Z,  t)  =  u0(X,  t)  +  Z<j>x{X,  t )  -  Z3Cl  (<j>x(X,  t)  +  ^  (la) 

w(X,Y,Z,t)  =w0(X,t)  (lb) 

where  the  X  coordinate  is  taken  along  the  beam  length,  the  Z  coordinate  along  the  thickness 

direction  of  the  beam,  Mo  is  the  axial  displacement  of  a  point  on  the  mid-plane  ( X ,  0,  0)  of  the 
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beam  and  Wq  represents  the  transverse  deflection  of  the  mid-plane.  When  the  deformation 
is  small  the  parameter  (j)x(X,t)  may  be  interpreted  as  the  rotation  of  the  transverse  normal. 
The  constant  c±  is  equal  to  C\  =  4/(3 h2),  where  h  is  the  height  of  the  beam  and  b  is  the 
beam  width.  The  displacement  field  of  the  Reddy  beam  theory  suggests  that  a  straight  line 
perpendicular  to  the  undeformed  mid-plane  becomes  a  cubic  curve  following  deformation,  as 
can  be  seen  in  Fig.  1. 


FIGURE  f.  Deformation  of  a  beam  structure  according  to  the  third-order 
Reddy  beam  theory:  (a)  undeformed  configuration  and  (b)  deformed  configu¬ 
ration. 


2.2.  The  effective  strain  tensor  for  the  simplified  theory.  In  the  mechanical  analysis 
of  deformable  solids,  it  is  necessary  to  employ  stress  and  strain  measures  that  are  consistent 
with  the  deformations  realized  [34,  6].  When  the  deformations  of  the  body  are  large,  there 
are  a  variety  of  strain  measures  that  may  be  employed.  In  our  formulation  we  employ  a  total 
Lagrangian  description  of  the  deformation.  In  such  analysis,  the  Green- Lagrange  strain 
tensor  E  constitutes  an  appropriate  measure  of  the  strain  at  a  point  in  the  body.  For  the 
present  analysis  the  non-zero  components  of  E  may  be  expressed  as 


E; xx  = 


du 


dX 


du\  f  dw 

dx)  +  VftX 

1  f  du  dw  du  du 


2i 


Exz  2  \az  +  ax  +  axaz 


(2a) 

(2b) 

(2c) 


2  \  ,/Z 

In  the  present  formulation  we  wish  to  develop  a  finite  element  framework  that  is  applicable 

under  loading  conditions  that  produce  large  transverse  displacements,  moderate  rotations 
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(10  —  15°)  and  small  strains  [32],  Under  such  conditions  it  is  possible  to  neglect  the  underlined 
terms  in  the  above  definition  of  the  Green- Lagrange  strain  tensor.  Consequently,  we  employ  a 
reduced  form  of  the  Green-Lagrange  strain  tensor,  denoted  by  £,  whose  non- zero  components 
may  be  expressed  as 


where  C2  =  3ci.  The  strain  components  associated  with  the  linearized  strain  tensor  e  are 
commonly  called  the  the  von  Karman  strain  components.  For  a  comparison  of  numerical 
results  obtained  using  the  above  simplified  theory  with  the  full  nonlinear  theory  for  clastic 
structures,  we  refer  to  the  work  of  Ba§ar  et  al.  [5].  It  is  important  to  note  that  the  material 
coordinates  appearing  in  the  definition  of  the  reduced  strain  components  and  throughout  the 
remainder  of  this  paper  are  denoted  as  ( x ,  y,  z )  as  a  reminder  that  the  present  formulation 
is  applicable  to  small  strains  and  moderate  rotations,  and  is  therefore  a  linearization  of  the 
more  general  finite  deformation  theory. 


2.3.  Linear  viscoelastic  constitutive  equations.  For  linear  viscoelastic  materials,  the 
constitutive  equations  relating  the  components  of  the  second  Piola  Kirchhoff  stress  tensor  S 
to  the  Green-Lagrange  strain  E  may  be  expressed  in  terms  of  the  following  set  of  integral 
equations 

S (t)  =  <C(0)  :  E (t)  +  [  C (t-s):  E (s)ds  (4) 

Jo 

where  C (t  —  s)  =  dC(t  —  s)/d(t  —  s )  and  C (t)  is  the  fourth-order  viscoelasticity  relaxation 
tensor.  In  the  present  analysis  the  above  expression  reduces  to 


&xx (x,  t)  =  E(0)£xx(x,t)  + 

<7xa(x,i)  =  G(0)7x*(x,t)  + 


E(t  -  s)£xx(x,s)ds 


G(t  -  s)'fxz(x,s)ds 


(5a) 

(5b) 


where  axx  and  axz  are  the  nonzero  components  of  second  Piola  Kirchhoff  stress  tensor  used  in 

the  present  simplified  formulation.  The  quantities  E(t)  and  G(t)  are  the  relaxation  moduli. 

The  specific  forms  of  E(t)  and  G(t)  will  depend  upon  the  material  model  employed.  For  the 

present  analysis  we  assume  that  these  relaxation  functions  can  be  expanded  as  Prony  series 
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of  order  NPS  as 


NPS  NPS 

E(t)  =  E0  +  ^2El(t),  G(t)  =  G0  +  J2Gi(t)  (6) 

1=1  1=1 

where  E[(t)  and  Gi(t)  have  been  defined  as  (following  the  generalized  Maxwell  model) 

Ei(t)  =  Eie-t/T",  Gi(t)  =  Gte-^F  (7) 

It  is  important  to  emphasize  that  the  Prony  series  representation  of  the  viscoelastic  relax¬ 
ation  moduli  is  critical  for  the  implementation  of  efficient  temporal  numerical  integration 
algorithms  of  the  integral  type  viscoelastic  constitutive  equations  considered  in  this  paper. 
We  note  in  passing  that  effective  temporal  integration  algorithms  for  alternative  classes 
of  viscoelastic  constitutive  equations,  such  as  fractional  derivative  models,  have  also  been 
adopted  in  the  literature  (see  for  example  Refs.  [10,  11,  12,  15,  42]). 


3.  Weak  form  Galerkin  finite  element  model 


3.1.  Weak  formulation.  The  weak  form  Galerkin  finite  element  model  of  the  third  order 
Reddy  beam  theory  may  be  developed  by  applying  the  principle  of  virtual  displacements  to 
a  typical  beam  as  viewed  in  the  reference  configuration.  The  dynamic  form  of  the  virtual 
work  statement  may  therefore  be  expressed  as 


5W  =  -SJC  +  6W  i  +  6We 


,dS 


=  f  (5u  •  p0ii  +  SE  :S-5u-  p0b)dV  -  f  hu  •  tc 
Jb  J  r  fj 

=  f  f  (<5u  •  p0ii  +  5e  :  cr  —  hu  •  p0h)dAdx  —  f  <5u  •  t0dS  =  0 


(8) 


where  6IC  is  the  virtual  kinetic  energy,  hWi  is  the  internal  virtual  work  and  JWe  is  the 
external  virtual  work.  The  additional  quantities  po,  b  and  to  are  the  density,  body  force  and 
traction  vector,  respectively.  The  above  expression  constitutes  the  weak  form  of  the  classical 
Euler-Lagrange  equations  of  motion  of  a  continuous  body. 

In  the  finite  element  method  we  assume  that  the  material  domain  =  [0,  L\  is  partitioned 
into  a  set  of  NE  non-overlapping  sub-domains  f le  =  called  finite  elements,  such  that 

O  =  U^=i  ^e-  The  resulting  variational  problem  associated  with  the  weak  formulation  of 
the  Reddy  beam  equations  may  therefore  be  expressed  as  follows:  find  ( uo,wo,(j)x )  G  V  such 


that  for  all  (5uq,  5wq,  5<px)  G  X  the  following  expressions  hold  within  each  element: 


0 

0 


0 


(c)Su  \ 

IqSuqUo  +  -^-Nxx  -  Suof  )  dx  -  5u0(xa)Q  1  -  5u0(xb)Q5 


rxb 


I0Sw0w0  + 


05  w, 


0  /  „2 


'  Xa  L 

d25w, 


Ox 


c\h 


dido 

Ox 


-  J 


4ycr 


+ 


05wo  f  Owq 


dx  V  Ox 


P.r.r  Qx 


Ox 2 


Ql 


CXb 


CiPxx  ~  5w0q 


05wr 


dx  -  Q25w0(xa)  -  Q65w0(xb )  -  Q3 


05  wc 
Ox 


Ox 


X=Xb 


Own 


Ox 


05(j)x 


54>x  - h  4>x  )  +  5<px{Qx  —  C2Rx)  H - ^ - ( Mxx  —  C\PXX ) 


Ox 


dx 


Q  idcpxiXa)  Q&dc^xiXh) 


(9a) 

(9b) 


(9c) 


where  V  and  X  are  appropriate  function  spaces.  It  is  important  to  note  that  in  the  finite  ele¬ 
ment  implementation,  wq  must  be  approximated  using  Hcrmite  polynomials.  Consequently, 
the  nodal  degrees  of  freedom  of  the  resulting  finite  element  model  will  contain  not  only  wq 
but  also  its  derivative.  For  the  sake  of  brevity  we  have  omitted  the  superscript  e  from  quan¬ 
tities  appearing  in  the  above  equations  (e.g.,  xa  and  xb)  and  throughout  the  remainder  of 
this  work.  The  quantities  /  and  q  appearing  above  are  the  distributed  axial  and  transverse 
loads  respectively.  We  have  also  introduced  the  following  constants 

Ii  =  PoDi  =  po  /  zldA ,  J4  =  —  c\Iq),  K2  =  I2  —  2C1/4  +  c\Iq  (10) 

J  A 

The  internal  stress  resultants  Nxx,  Mxx,  Pxx,  Qx  and  Rx  are  defined  as 


and  can  be  expressed  in  terms  of  the  generalized  displacements  (uq,  wq,  4>x)  through  the  use 
of  the  viscoelastic  constitutive  equations.  The  quantities  Nxx,  Mxx  and  Qx  are  the  internal 
axial  force,  bending  moment  and  shear  force.  In  addition,  Pxx  and  Rx  arc  higher  order  stress 
resultants  that  arise  in  the  third  order  beam  theory  due  to  the  cubic  expansion  of  the  axial 
displacement  field.  The  quantities  Q3  (j=l,..,8)  are  the  externally  applied  generalized  nodal 
forces. 
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3.2.  Semi-discrete  finite  element  equations.  In  this  subsection  we  develop  the  semi¬ 
discrete  finite  element  equations  associated  with  the  third  order  Reddy  beam  theory.  Within 
a  typical  finite  element  the  dependent  variables  may  be  adequately  approximated  using  the 
following  interpolation  formulas 


u0  (x,t)  =  ^2 


2  n 


3= 1 


W0{x,t)  =  ^2  A f{t)^f\x) 

3= 1 


(12) 


&,t)  =  J2 


X) 


3= 1 


where  a  space-time  decoupled  formulation  has  been  adopted  and  n  represents  the  number 
of  nodes  per  element.  In  the  finite  element  method,  the  geometry  of  each  element  is  char¬ 
acterized  using  the  standard  isoparametric  mapping  from  the  master  element  Vte  =  [—1,  +1] 
to  the  physical  element  0e  =  [x®,  x|],  The  quantities  A f\t),  A f\t)  and  A f\t)  are  the 
generalized  displacements  at  the  element  nodes.  In  addition  3j/p(x)  and  x )  are  the 

(n  —  l)th-order  Lagrange  and  (2 n  —  l)th-order  Hermite  interpolation  functions  respectively. 
Inserting  Equation  Eq.  (12)  into  Eqs.  (9a)  through  (9c)  results  in  the  semi-discrete  finite 
element  equations  for  the  RBT.  The  resulting  set  of  equations  for  a  typical  element  may  be 
expressed  as  at  the  current  time  t  as 


[ilT]{Ae}  +  [Ke}{Ae}  +  J  {A e(t,s)}ds  =  {Fe} 


(13) 


The  element  level  equations  may  be  partitioned  into  the  following  equivalent  set  of  expres¬ 
sions  ^ 

[MaP]{A^}  +  [A'q/3]{A(/3)}  +  [  {A{a\t,s)}ds  =  {F(a)}  (14) 

Jo 

where  a  and  f3  range  from  1  to  3  and  Einstein’s  summation  convention  is  implied  over  f3. 
Expressions  for  determining  the  components  of  the  partitioned  element  coefficient  matrices 
and  vectors  are  given  explicitly  in  Appendix  A. 


3.3.  Fully-discrete  finite  element  equations.  In  this  subsection  we  develop  the  fully 
discretized  finite  element  equations  for  the  Reddy  beam  theory.  We  begin  by  partitioning 
the  time  interval  [0,r]  C  K  of  interest  in  the  analysis,  where  r  >  0,  into  a  set  of  N  non- 
overlapping  subintervals  such  that  [0,  r]  =  (J^1Xfc,  and  Ik  =  [tk,tk+ 1]-  The  solution  may 
then  be  obtained  incrementally  by  solving  an  initial  value  problem  within  each  subinter¬ 
val  Zfc,  where  we  assume  that  the  solution  is  known  at  t  —  tk ■  Within  each  subregion  it  is 


therefore  necessary  to  introduce  approximations  for  both  the  temporal  derivatives  of  the  gen¬ 
eralized  displacements  (resulting  from  the  inertia  terms)  as  well  as  the  convolution  integrals 
(resulting  from  the  viscoelastic  constitutive  model  of  the  material).  Since  temporal  integra¬ 
tion  of  the  inertia  terms  is  relatively  straightforward,  we  restrict  the  current  discussion  to 
discretization  of  the  quasi-static  form  of  the  semi-discrete  finite  element  equations  only.  In 
this  work,  we  approximate  the  convolution  integrals  using  the  trapezoidal  rule  within  each 
time  subinterval.  We  further  introduce  a  two-point  recurrence  formula,  associated  with  a  set 
of  history  variables  that  are  evaluated  at  the  quadrature  points,  that  is  utilized  to  effectively 
advance  the  numerical  solution  from  one  time  step  to  the  next  such  that  data  history  is 
only  necessary  from  the  immediate  previous  time  step.  Although  not  entirely  the  same,  the 
adopted  procedure  has  its  roots  in  many  of  the  key  ideas  presented  in  the  pioneering  work 
of  Taylor  et  al.  [38] ,  where  the  finite  element  method  was  first  employed  to  solve  problems 
in  viscoelasticity  using  algorithms  based  on  recurrence  formulas. 

We  assume,  without  loss  of  generality,  that  the  quasi-static  semi-discrete  finite  element 
equations  have  been  successfully  integrated  temporally  up  until  t  =  tk.  Our  goal,  therefore, 
is  to  numerically  integrate  the  finite  element  equations  over  the  subinterval  Xk  to  obtain  the 
solution  for  the  generalized  displacements  at  t  —  tk+ 1-  Before  proceeding  we  must  empha¬ 
size  that  all  subsequent  discussions  regarding  efficient  recurrence  based  temporal  integration 
strategies  rely  on  the  following  multiplicative  decompositions  of  the  Prony  series  terms  ap¬ 
pearing  in  the  definition  of  the  relaxation  moduli  [37] 

4(4+ 1  -  s)  =  e-Atfc+i/^(4  -  ^),  4(4+1  -  s)  =  e-A^A ? 6t(tk  -  s )  (15) 

where  A tk+i  =  4+i  —  4  is  the  time  step  associated  with  subinterval  Xk.  With  the  above 
formulas  in  mind,  we  note  that  the  components  of  A“(4+i,  s)  may  be  conveniently  expressed 
as 

XI OL 

A"(4+i,s)  =  ^JA  f(4+1,s)  (16) 

3=1 

where  ri\  =  1,  ri2  =  3  and  713  =  2.  The  components  -^“(4+1,5)  can  be  decomposed 
multiplicatively  using  the  following  general  formula 

NGP  NPS 

JA?(  4+,,s)  =  (17) 

m=  1  /=! 
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In  the  above  expression  we  have  employed  the  Gauss-Legendre  quadrature  rule  in  evaluation 
of  all  spatial  integrals  (resulting  in  summation  over  m).  The  quantity  Wm  represents  the 
mth  quadrature  weight  associated  with  the  Gauss-Legendre  quadrature  rule.  Summation 
over  l  is  due  to  the  Prony  series  representation  of  the  relaxation  moduli.  The  multiplicative 
decomposition  of  each  3  A “(ifc+1,s)  is  essential  for  the  recurrence  based  integration  strategy. 
The  components  of  Jmxf(tk+ 1)  are  used  to  store  the  discrete  finite  element  test  functions  as 
well  as  any  nonlinear  quantities  associated  with  the  first  variation  of  the  simplified  Green- 
Lagrange  strain  tensor.  In  the  present  formulation  the  components  of  3mxf{tk+i)  are  defined 
as 


mXi  (^fe+l) 

mXi  (^fe+l) 

2mXKtk+ 1) 

mXlih+l) 

mXi(tk+ 1) 


d^1\xm) 


dx 


dw0(xm,tk+ 1)  d^f\xr 


dx 

d2ll)f\xm) 

dx2 


dx 


d^f\xm) 

dx 

rrtX.i 


mXUtk+l)  =  ♦f’W 


(18a) 

(18b) 

(18c) 

(ISd) 
( 18e) 

(18f) 


In  the  above  expression,  xm  represents  the  value  of  x  as  evaluated  at  the  mth  quadrature 
point  of  a  given  finite  element.  The  isoparametric  mapping  Vte  fle  used  to  characterize  the 
geometry  of  each  element  allows  for  simple  evaluation  of  such  expressions.  The  components 
of  j/3a(Atk+i)  are  defined  as 

}^{Atk+ 1)  =  ]/32(Atk+i)  =  2f32(Atk+1)  =  }/33(Atk+i)  =  e-Atk+l/TF 


(19) 

f/3\Atk+1)  =  2f33(Atk+1)  =  e~Atk+1^T° 

Likewise,  the  components  of  jmKa(tk,s)  may  be  determined  using  the  following  formulas 

(20a) 
(20b) 
(20c) 


imKl(*fc>s)  =  El(fk  ~  S)D0 

ImA  {tki  s)  =  lmK  {tki  s ) 

K2(tk,s)  =  Ei(tk  -s)(  c(D6 


du0(xm,s )  1/ 

^  dvJo(xm,  s')  \  2 

dx  2 ' 

\  dx  ) 

2  ..2 
Im 


2r,  d2w0(  s  )  r  9<j)x(  '-0 


dx 2 


~  La 


dx 


to 


lm,K  (tk,  s)  —  Gl(tk 
Im K  (tk,s)  =  El(tk 


s)As 


f  dw0(xm,s ) 
\  dx 


s) 


0(f) x  (xmi  s) 
Ox 


T  d2w0(xm,s)\ 

4  Ox2  J 


L«3(*fc,s)  =  fmK2(tk,s) 


(20d) 

(20e) 

(20f) 


It  is  important  to  note  that  the  components  of  3lmKa(tk,s)  have  been  defined  such  that  the 
following  multiplicative  recurrence  formulas  holds 


ImK“fe+ 1.S)  =  ^“(A4+l 


(21) 


The  above  expressions  are  admissible  on  account  of  the  assumption  that  the  relaxation 
parameters  are  expressed  in  terms  of  Prony  series. 

We  assume  that  at  t  =  tk  the  components  of  the  following  expression  are  known 

ntu  ngp  NPS 


/  jK(tk,s)ds  =  J2J2L 


X, 


'ML*' 


( tk)W„ 


(22) 


o 


m=  1  1=1 


where  3lmXa(tk)  is  a  set  of  history  variables  (stored  at  the  quadrature  points  of  each  element) 
that  are  of  the  form 

L  Xa(tk)=  [ tkimKa(tk,s)ds  (23) 

Jo 

We  note  that  jmXa( 0)  =0.  At  t  =  the  above  history  variables  are  known  and  there  is  no 
need  to  explicitly  evaluate  the  expression  appearing  on  the  right  hand  side  of  Eq.  (23).  At 
the  subsequent  time  step  t  =  tfc+1  Eq.  (22)  may  be  written  as 


fifc+l 


rtk 


ftfc+i 


3Ai(tk+1,s)ds  =  /  3Af(tk+i,s)ds+  /  3A?(tk+i,s)ds 

Jo  J  tk 

NGP  NPS 

=  ££'»■  x?(4+i);/9“(A4+i);mx“(4)iG. 

m.=l  Z=1 
rtk+l 

+  /  JA?(ifc+1,s)ds 

Jtk 


(24) 


It  is  important  to  note  that  we  have  expressed  the  first  integral  on  the  right  hand  side  of 
the  above  equation  in  terms  of  jmXa(tk )  (which  is  known  from  the  previous  time  step). 
To  integrate  the  remaining  expression  in  Eq.  (24)  over  the  subinterval  Xk  we  employ  the 


ll 


trapezoidal  rule  which  may  be  expressed  as 


fifc+l 


At 


3Ai(tk+l,s)ds  =  '^+1  pA“(4+i,4)  +  3Af(tk+i,tk+i)_ 


'tk 


At 


NGP  NPS 


k+ 1 


E5A  xt(tk+1y^a{Atk+1)  \jmKa(tk,tk ) 


m= 1  /=1 

+/tok<*(4,  4+i)]bhm 

As  a  result,  Eq.  (24)  can  be  written  in  the  following  simplified  form 

tk+1  NGP  NPS 

/  jK(tk+l,s)ds  =  X) 

*'°  m=l  i=l 


(26) 


where 


Zra 


fc+u 


_  Atk+ 1  j 


/3"(A4+1)  4)  +  imKa(tk,  tk+ 1) 


+  ^“(A4+i)L^a(4) 


(27) 


As  a  result,  in  Eq.  (26)  we  have  developed  a  general  expression  for  integrating  the  viscoelastic 
terms  up  to  any  discrete  instance  in  time.  The  expression  relies  on  a  recurrence  relationship 
defined  in  terms  of  the  set  of  history  variables  JlrnXa(tk+i).  These  variables  must  be  stored 
in  memory  at  the  immediate  previous  time  step  and  may  be  updated  to  the  subsequent 
time  step  in  accordance  with  the  procedure  outlined  in  Eq.  (27).  The  history  variables  are 
expressed  explicitly  in  Appendix  B.  It  is  now  possible  to  express  the  fully  discretized  finite 
element  equations  at  the  current  time  step  as 

[Ke]k+1{Ae}k+1  =  {Fe}k+1  -  {&}k+1  (28) 

where  each  of  the  components  in  the  above  equation  are  given  in  Appendix  B. 


3.4.  Iterative  solution  of  nonlinear  equations  using  Newton’s  method.  The  fully 

discretized  finite  element  equations  are  nonlinear  due  to  the  use  of  the  von  Karman  strain 
components  in  the  definition  of  the  effective  strain  tensor  e.  In  our  work,  we  adopt  the 
Newton  procedure  in  the  iterative  solution  of  the  nonlinear  finite  element  equations.  The 
resulting  linearized  finite  element  equations  are  of  the  form 


rthuAxy  =  -([A-itMA'Ci  -  {**}& + mho 
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i  (r+1) 


lb) 


ib) 


ib) 


ib) 


(29) 


where  represents  the  incremental  solution  at  the  (r  +  l)th  nonlinear  iteration. 

The  total  global  solution  at  the  (r  +  l)th  iteration  is  obtained  as 


{A}'-!1’  =  {iA}<;+1I)  +  {A}w 


k+ 1 


(30) 


The  element  tangent  stiffness  coefficient  matrix  [Te]^\  appearing  in  the  Newton  linearization 
of  the  finite  element  equations  may  be  expressed  (using  Einstein’s  summation  convention  over 
n)  as 

dQt 


d  Ke 

Te  =  Ke.  +  — — Ae  + 

13  13  dAe  n  dAe 


(31) 

All  quantities  comprising  the  tangent  stiffness  matrix  are  formulated  using  the  solution  from 
the  rth  iteration.  The  partial  derivatives  are  taken  with  respect  to  the  solution  at  the  current 
time  step.  Formulas  for  evaluating  the  components  of  the  element  tangent  stiffness  coefficient 
matrix  are  given  in  Appendix  C. 


3.5.  Finite  element  interpolation  of  dependent  variables.  It  is  well-known  that  low- 
order  finite  elements  for  beams  are  prone  to  locking  [34,  31,  33]  when  quadrature  rules 
are  employed  that  result  in  exact  integration  of  the  element  coefficient  matrices  and  force 
vectors.  To  circumvent  the  locking  phenomena,  we  consider  two  philosophically  dissimilar 
numerical  procedures.  In  the  first  approach,  we  employ  the  lowest  order  element  admissible 
in  the  formulation  (i.e.,  a  two-noded  element).  Selective  full  and  one  point  Gauss- Legendre 
quadrature  rules  are  applied;  where  reduced  integration  techniques  are  employed  on  all 
nonlinear  expressions  associated  with  the  finite  element  model.  This  element  is  denoted  as 
an  RBT-2-R  element  (meaning  a  two-noded  reduced  integration  RBT  element).  It  is  worth 
noting  that  this  element  requires  a  splitting  of  the  history  variables  into  subsets  associated 
with  the  full  and  reduced  quadrature  points.  In  the  second  approach,  we  construct  the  Reddy 
beam  finite  elements  using  high-polynomial  order  expansions  of  the  dependent  variables,  by 
systematically  increasing  the  number  of  nodes  per  finite  element.  In  this  approach,  the  same 
quadrature  formulas  may  be  employed  in  the  evaluation  of  all  expressions  appearing  in  the 
coefficient  matrices  and  force  vectors  of  the  finite  element  model.  The  resulting  elements 
are  denoted  in  this  work  as  RBT-n  elements,  where  n  represents  the  number  of  nodes  per 
element. 

In  the  proposed  high-order  finite  element  formulation,  we  employ  an  unequal  spacing  of 

the  nodes  within  each  element.  We  define  the  nodal  points,  within  the  master  element 
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Cte  =  [—1,  +1],  in  terms  of  the  roots  of  the  following  expression 


K  -  rn  +  !)£;«)  =  0  in  n' 


(32) 


where  Lp  (£)  is  the  Legendre  polynomial  of  order  p  [17]  and  n  —  p+1  represents  the  number 
of  nodes  per  element.  As  a  result,  the  quantities  where  i  —  1,  ..,n,  are  the  nodal  points 
associated  with  f 2e.  In  our  formulation  we  construct  the  high-order  Lagrange  interpolation 
functions  in  accordance  with  the  following  expression 


,,,  K-DK  +  i  )L'„(Q 

Wi  K>  p(p+i)Lr(m-(.) 


in 


(33) 


The  above  interpolation  functions  are  often  called  spectral  nodal  interpolation  functions  in 
the  literature  [20].  The  high-order  Hermite  interpolation  functions  ^•2')  associated  with  the 
master  element  fle  may  be  calculated  as 


</?’«)= Ter1?-1 


2  n 


3= 1 


(34) 


The  coefficients  d  1  appearing  in  the  above  expression  may  be  determined  by  imposing  the 


following  compatibility  conditions  on  the  interpolation  functions 


V4?-i  (0)  = 


#2? 


=  5. 


1 


IJ, 


d( 


=  d?(0)  =  0 


(35) 


where  i  and  j  both  range  from  1  to  n.  The  Hermite  interpolation  functions  ?/y  ;  associated 
with  the  physical  element  h2e  may  be  determined  as 


d7,K)  =  d7,K).  d?K)  =  ^d?K) 


(36) 


where  Je  is  the  Jacobian  of  the  element  coordinate  transformation  Vte  f2e.  The  above  for¬ 
mulas  may  be  utilized  to  generate  high-order  Lagrange  and  Hermite  interpolation  functions 
for  any  number  of  nodes  per  element.  The  standard  lowest  order  two-noded  element  may 
be  obtained  as  a  special  case.  The  interpolation  functions  associated  with  a  six-node  RBT 
finite  element  are  shown  in  Fig.  2. 


4.  Numerical  examples 

In  this  section,  numerical  results  are  presented  and  tabulated  for  the  mechanical  response 

of  viscoelastic  beam  structures  obtained  using  the  proposed  finite  element  formulation  for 
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FIGURE  2.  Interpolation  functions  for  a  high-order  RBT  finite  element  where 
n  =  6  and  i  =  1  (a)  Lagrange  interpolation  functions  (b)  Hermite 

interpolation  functions  fyi-i  and  (c)  Hermite  interpolation  functions  ■ 


the  third-order  Reddy  beam  theory.  The  results  presented  in  this  section  have  been  obtained 
using  the  Newton  solution  procedure  described  in  the  previous  section.  Nonlinear  conver¬ 
gence  is  declared  at  time  £*.+i  once  the  Euclidean  norm  of  the  normalized  difference  between 
the  nonlinear  iterative  solution  increments  (i.e. ,  [[{A}^^  —  {A}M||/||{A}<yi1,||),  is  less 
than  10-6. 

The  material  model  utilized  in  the  quasi-static  numerical  studies  is  based  upon  the  exper¬ 
imental  results  tabulated  by  Lai  and  Bakker  [23]  for  a  glassy  amorphous  polymer  material 
(PMMA).  The  Prony  series  parameters  for  the  viscoelastic  relaxation  modulus  given  in  Table 
1  were  calculated  by  Payette  and  Reddy  [29]  from  the  published  compliance  parameters  [23] . 

As  in  the  work  of  Chen  [8]  and  Payette  and  Reddy  [29],  we  assume  that  Poisson’s  ratio  is 
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(37) 


time  invariant.  As  a  result,  the  shear  relaxation  moduli  is  given  as 


G{t) 


m 

2(l  +  i/) 


where  Poisson’s  ratio  is  taken  to  be  v  =  0.40  [22], 


Table  1.  Viscoelastic  relaxation  parameters  for  a  PMMA. 


E0 

205.7818  ksi 

E\ 

43.1773  ksi 

rf 

9.1955  x  10"1  sec. 

E2 

9.2291  ksi 

t# 

9.8120  x  10°  sec. 

e3 

22.9546  ksi 

te 

'3 

9.5268  x  101  sec. 

e4 

26.2647  ksi 

te 

'4 

9.4318  x  102  sec. 

e5 

34.6298  ksi 

te 

1 5 

9.2066  x  103  sec. 

Eq 

40.3221  ksi 

te 

'6 

8.9974  x  104  sec. 

e7 

47.5275  ksi 

te 

'7 

8.6852  x  105  sec. 

E8 

46.8108  ksi 

te 

'8 

8.5143  x  106  sec. 

Eg 

58.6945  ksi 

te 

'9 

7.7396  x  107  sec. 

4.1.  Quasi-static  deflection  of  various  thin  beams  under  uniform  loading.  In  this 
first  example  we  consider  a  viscoelastic  beam  of  length  L  =  100  in.  and  cross  section  1  in.  x 
1  in.  At  t  —  0  sec.  the  beam  is  subjected  to  a  uniform  vertically  distributed  load  go  =  0.25 
lbf/in.  Due  to  symmetry  about  x  =  L/ 2,  it  is  only  necessary  to  model  half  of  the  physical 
domain.  To  assess  the  performance  of  various  finite  element  discretizations  in  circumventing 
the  locking  phenomena,  we  consider  the  following  three  sets  of  boundary  conditions: 


(1)  Hinged  at  both  ends 

m(0,t)  =  u0(L/2,t)  =  -^-(L/2,t)  =  <f>x(L/2,t)  =  0 

(2)  Pinned  at  both  ends 

uo{0,t)  =  wo(0,t)  =  u0(L/2,t)  =  -^-(L/2,t)  =  c i>x(L/2,t )  =  0 

(3)  Clamped  at  both  ends 

■u0(0,  t )  =  w0(0,  t)  =  -7^(0,  t)  =  ^(0,  t)  =  0 

r\ 

ML/ 2,  t)  =  ~^r(L/  2,  t)  =  ML/ 2,  t)  =  o 


(38) 


(39) 


(40) 
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(b) 


time,  t 


Figure  3.  Maximum  vertical  deflection  wo(L/2,t)  of  various  viscoelastic 
beams,  each  subjected  to  a  uniform  vertically  distributed  load  go-  2  RBT- 
6  elements  employed  in  each  finite  element  discretization:  (a)  hinged-hinged 
beam  configuration  and  (b)  pinned-pinned  and  clamped-clamped  beam  con¬ 
figurations. 


17 


TABLE  2.  Comparison  of  the  quasi-static  finite  element  solutions  for  the  maxi¬ 
mum  vertical  deflection  w0(L/ 2,  t )  of  viscoelastic  beams  subjected  to  a  uniform 
load  go  for  three  different  boundary  conditions. 


Time,  t 

RBT-2 

RBT-2-R 

RBT-3 

RBT-4 

RBT-6 

Hinged-hinged 

0 

5.4740 

7.2840 

7.2277 

7.2946 

7.2980 

200 

6.1234 

8.6052 

8.5170 

8.6169 

8.6217 

400 

6.1895 

8.7473 

8.6552 

8.7592 

8.7641 

600 

6.2295 

8.8340 

8.7396 

8.8460 

8.8510 

800 

6.2615 

8.9035 

8.8071 

8.9156 

8.9207 

1,000 

6.2886 

8.9627 

8.8646 

8.9748 

8.9799 

1,200 

6.3119 

9.0138 

8.9143 

9.0259 

9.0311 

1,400 

6.3322 

9.0584 

8.9575 

9.0705 

9.0758 

1,600 

6.3499 

9.0975 

8.9956 

9.1098 

9.1150 

1,800 

6.3656 

9.1322 

9.0293 

9.1445 

9.1498 

Pinned-pinned 

0 

1.2442 

1.2493 

1.2452 

1.2452 

1.2452 

200 

1.3233 

1.3291 

1.3242 

1.3242 

1.3242 

400 

1.3313 

1.3371 

1.3322 

1.3322 

1.3322 

600 

1.3361 

1.3420 

1.3370 

1.3370 

1.3370 

800 

1.3399 

1.3459 

1.3409 

1.3409 

1.3409 

1,000 

1.3432 

1.3492 

1.3441 

1.3441 

1.3441 

1,200 

1.3460 

1.3520 

1.3470 

1.3470 

1.3469 

1,400 

1.3484 

1.3545 

1.3494 

1.3494 

1.3494 

1,600 

1.3506 

1.3566 

1.3515 

1.3515 

1.3515 

1,800  1.3525 

Clamped- clamped 

1.3585 

1.3534 

1.3534 

1.3534 

0 

0.9037 

0.9098 

0.9106 

0.9108 

0.9109 

200 

0.9918 

0.9992 

0.9993 

0.9995 

0.9997 

400 

1.0007 

1.0082 

1.0083 

1.0084 

1.0086 

600 

1.0060 

1.0136 

1.0136 

1.0138 

1.0140 

800 

1.0103 

1.0180 

1.0180 

1.0182 

1.0183 

1,000 

1.0139 

1.0216 

1.0216 

1.0218 

1.0220 

1,200 

1.0170 

1.0248 

1.0247 

1.0249 

1.0251 

1,400 

1.0197 

1.0275 

1.0275 

1.0277 

1.0278 

1,600 

1.0221 

1.0299 

1.0298 

1.0300 

1.0302 

1,800 

1.0242 

1.0321 

1.0319 

1.0322 

1.0323 

In  Table  2  and  Fig.  3  we  summarize  numerical  results  for  the  maximum  vertical  deflection 

of  the  viscoelastic  beam  for  the  three  different  sets  of  boundary  conditions  listed  above.  The 

tabulated  results  were  obtained  using  10  RBT-2  elements  (11  nodes),  10  RBT-2-R  elements 

(11  nodes),  5  RBT-3  elements  (11  nodes),  3  RBT-4  elements  (10  nodes)  and  2  RBT-6 
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elements  (11  nodes).  An  equal  time  increment  of  At  =  1.0  sec.  was  employed  for  all  time 
steps.  To  insure  convergence  of  the  nonlinear  solution  procedure,  the  instantaneous  clastic 
solution  (at  1  =  0  sec.)  was  obtained  through  the  use  of  five  load  steps.  At  each  subsequent 
time  step,  the  finite  element  equations  were  solved  iteratively  using  the  Newton  procedure 
without  employing  load  stepping.  Typically,  only  2  or  3  Newton  iterations  were  necessary 
to  satisfy  the  nonlinear  convergence  criterion. 

It  is  interesting  to  note  that  the  numerical  results  for  all  finite  element  discretizations  are 
comparable  with  the  exception  of  the  case  where  the  beam  is  subjected  to  hinged  boundary 
conditions  at  both  ends.  For  this  problem,  the  RBT-2  finite  element  clearly  suffers  from 
locking.  It  is  evident,  however,  that  polynomial  refinement  of  the  solution  naturally  alleviates 
the  locking.  In  fact,  the  RBT-3  element  is  almost  completely  locking  free.  Overall  we  find 
that  this  element  is  less  prone  to  locking  (when  n  is  small)  than  the  Timoshenko  beam 
element  employed  previously  by  Payette  and  Reddy  [29].  The  results  obtained  for  the  RBT- 
6  element  are  spatially  fully  converged  and  compare  well  with  the  Timoshenko  beam  results 
obtained  using  high-order  polynomial  expansions  [29]. 

For  the  hinged-hinged  beam  configuration,  the  vertical  deflection  coincides  with  the  exact 
solution  of  the  geometrically  linear  theory.  In  Table  3  we  compare  numerical  results  obtained 
using  2  RBT-6  beam  elements  with  the  exact  solution  for  the  Timoshenko  beam  theory  given 
by  Fliigge  [14]  as 


w0(L/2,t)  = 


5  q0L4 
384 Do 


r, ,  8(i +") 

5  K 

w  J 

D(t) 


(41) 


where  D{t)  is  the  creep  compliance  and  n  is  the  shear  correction  factor.  The  error  in  the 
numerical  solution  due  to  temporal  integration  via  the  trapezoidal  rule  tends  to  over-predict 
the  deflection  of  the  beam  as  is  evident  in  Table  3  (where  numerical  solutions  obtained  using 
various  time  increment  sizes  are  compared). 


4.2.  Quasi-static  deflection  of  a  thick  beam  under  uniform  loading.  In  this  next 

example  we  consider  a  thick  (i.e.,  L/h  <  20)  viscoelastic  beam  to  demonstrate  the  ability  of 

the  Reddy  beam  finite  element  formulation  to  accurately  account  for  deformations  associated 

with  shearing.  We  modify  the  thin  beam  problem  of  the  previous  example  by  letting  L  =  10 

in.,  q  =  25.0  lbf/in  and  At  =  1.0  sec.  All  other  geometric  and  material  parameters  are 

the  same  as  in  the  previous  example.  In  Table  4  numerical  results  are  presented  for  the 

transverse  deflection  of  pinned-pinned  and  clamped-clamped  beams.  The  same  number  of 

elements  (per  element  type)  are  employed  as  in  the  previous  example.  In  Table  4  we  also 
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Table  3.  Analytical  and  finite  element  solutions  for  the  maximum  quasi¬ 
static  vertical  deflection  wo(L/2,t)  of  a  hinged- hinged  beam  under  uniform 
transverse  loading  q0. 


Time,  t 

Maximum  vertical  deflection,  wo(L/2,t) 

Exact 

At  =  0.1 

At  =  1.0 

At  =  2.0 

At  =  5.0 

At  =  10.0 

0 

7.2980 

7.2980 

7.2980 

7.2980 

7.2980 

7.2980 

200 

8.5429 

8.5437 

8.6217 

8.8493 

10.2278 

14.7260 

400 

8.6827 

8.6835 

8.7641 

8.9993 

10.4291 

15.1493 

600 

8.7680 

8.7689 

8.8510 

9.0910 

10.5524 

15.4107 

800 

8.8364 

8.8372 

8.9207 

9.1645 

10.6513 

15.6214 

1,000 

8.8945 

8.8954 

8.9799 

9.2270 

10.7356 

15.8022 

1,200 

8.9448 

8.9456 

9.0311 

9.2811 

10.8087 

15.9597 

1,400 

8.9886 

8.9895 

9.0758 

9.3282 

10.8726 

16.0983 

1,600 

9.0271 

9.0280 

9.1150 

9.3697 

10.9288 

16.2210 

1,800 

9.0612 

9.0621 

9.1498 

9.4064 

10.9787 

16.3306 

compare  results  from  the  Reddy  beam  theory  with  finite  element  results  obtained  using  a 
low  order  reduced  integration  finite  element  model  based  on  the  Euler-Bernoulli  beam  theory 
(which  does  not  account  for  deformations  associated  with  shearing). 

Table  4.  Comparison  of  the  quasi-static  finite  element  solutions  for  the 
maximum  vertical  deflection  wo(L/2,t)  of  thick  pinned-pinned  and  clamped- 
clamped  viscoelastic  beams  under  uniform  transverse  loading  q0. 


Time  t 

Maximum  vertical  deflection,  wo(L/2,t) 

pinned-pinned 

clamped-clamped 

EBT-2-R 

RBT-2-R 

RBT-6 

EBT-2-R 

RBT-2-R 

RBT-6 

0 

0.07184 

0.07362 

0.07367 

0.01459 

0.01645 

0.01653 

200 

0.08437 

0.08643 

0.08649 

0.01724 

0.01943 

0.01952 

400 

0.08571 

0.08779 

0.08785 

0.01752 

0.01975 

0.01985 

600 

0.08652 

0.08862 

0.08869 

0.01769 

0.01995 

0.02004 

800 

0.08717 

0.08929 

0.08935 

0.01783 

0.02010 

0.02020 

1,000 

0.08773 

0.08985 

0.08992 

0.01795 

0.02024 

0.02033 

1,200 

0.08821 

0.09034 

0.09041 

0.01805 

0.02035 

0.02045 

1,400 

0.08862 

0.09076 

0.09083 

0.01814 

0.02045 

0.02055 

1,600 

0.08899 

0.09114 

0.09121 

0.01822 

0.02054 

0.02064 

1,800 

0.08931 

0.09147 

0.09154 

0.01829 

0.02062 

0.02072 

4.3.  Quasi-static  deflection  of  a  thin  beam  due  to  time-dependent  loading.  For  this 

next  example  we  employ  the  geometric  parameters,  material  properties  and  hinged-hinged 
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boundary  conditions  utilized  in  the  first  numerical  example.  We  replace  the  stationary 
uniformly  distributed  load  with  the  following  quasi-static  transverse  load 

q(t)  =  q0  | H(t)  -  ^  [(t  -  ar)H(t  -  ar )  -  (t  -  Pr)H(t  -  pr)\ |  (42) 

where  go  =  0.25  lbf/in,  r  =  200  sec.  and  H(t )  is  the  Heaviside  function.  The  parameters 
0  <  a  <  /3  <  1  are  constants  whose  values  may  be  appropriately  adjusted.  The  load  function 
above  is  constant  for  0  <  t  <  ar  and  decays  linearly  from  t  =  ar  to  t  =  (3t,  after  which  the 
load  is  maintained  at  zero.  We  utilize  the  above  loading  function  to  numerically  demonstrate 
that  the  finite  element  model  correctly  predicts  that  the  viscoelastic  beam  will  eventually 
recover  its  original  configuration  upon  removal  of  all  externally  applied  mechanical  loads. 
The  numerical  solution  for  the  problem  is  presented  in  Figure  4  for  various  values  of  a  and 
p.  It  is  clear  that  in  all  cases,  the  beam  tends  to  recover  its  original  configuration  as  t  tends 
to  infinity. 


Figure  4.  Maximum  vertical  deflection  wo(L/2,t)  of  a  hinged-hinged  vis¬ 
coelastic  beam  subjected  to  a  time-dependent  transverse  load  q(t). 

4.4.  Fully  dynamic  response  of  hinged  viscoelastic  beams.  As  a  final  example,  we 

consider  the  fully  transient  mechanical  response  of  viscoelastic  beams  under  mechanical 

loading.  For  this  example  we  employ  a  simple  three  parameter  solid  model  utilized  previous 
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by  Chen  [8].  In  the  standard  three  parameter  solid  model,  the  relaxation  modulus  may  be 
expressed  as 

E(t)  =  (l  -  e-^)  +  kie-'t'f  (43) 

k  1  +  rC2  V  ‘ 

where  in  the  present  example  k\  =  9.8  x  107  N/m2  and  k?  =  2.45  x  107  N/m2.  The  relaxation 
time  is  of  the  form  rf  =  rj/(ki  +  and  the  material  density  is  taken  to  be  p0  —  500  kg/m3. 
A  constant  Poisson  ratio  of  u  —  0.3  is  assumed. 

We  consider  a  beam  with  hinged  boundary  conditions  at  both  ends.  The  beam  length, 
width  and  thickness  are  taken  to  be  L  =  10  m,  b  —  2  m  and  h  =  0.5  m  respectively.  We 
consider  two  loading  scenarios.  In  loading  scenario  (1)  a  uniformly  distributed  transverse 
load  is  specified  along  the  entire  length  of  the  beam  as  q(t)  =  qoH(t )  N/m,  where  qo  =  10. 
Likewise,  for  loading  scenario  (2)  a  periodic  concentrated  force  is  applied  at  the  center  of 
the  beam  as  F(t)  =  qosm(nt)  N,  where  qo  =  50.  In  the  finite  element  discretization  of 
both  problems,  we  employ  2  RBT-6  elements  of  equal  size.  As  in  the  previous  examples, 
symmetry  is  once  again  exploited  in  construction  of  the  finite  element  meshes.  We  utilize  the 
Newmark-/?  procedure  [26]  for  performing  temporal  integration  of  the  inertia  terms  appearing 
in  the  fully  transient  beam  finite  element  equations.  The  Newmark  integration  parameters 
are  chosen  in  accordance  with  the  constant-average  acceleration  method  [35] .  Both  transient 
problems  are  solved  over  a  total  time  interval  of  20  sec.  For  loading  scenario  (1)  we  employ 
500  time  steps,  while  1,000  time  steps  are  utilized  for  loading  scenario  (2). 

Numerical  results  for  loading  scenarios  (1)  and  (2)  are  presented  in  Figs.  5  and  6  respec¬ 
tively.  In  each  figure  we  present  both  viscoelastic  and  clastic  solutions  (where  the  Young’s 
modulus  is  obtained  by  taking  FJelastic  =  A(0)).  As  expected,  the  viscoelastic  effects  tend  to 
add  damping  to  what  would  otherwise  be  purely  elastic  response.  In  Fig.  5  we  present  fully 
transient  viscoelastic  results  using  two  different  values  for  77.  This  problem  is  also  solved  using 
the  quasi-static  visoelasticity  solution  procedure.  It  is  evident  that  the  transient  viscoelas¬ 
tic  solution  approaches  the  steady-state  quasi-static  viscoelastic  response  once  a  sufficiently 
long  enough  period  of  time  has  transpired.  For  loading  scenario  (2)  the  viscoelasticity  effects 
clearly  reduce  the  overall  amplitude  of  the  forced  vibrational  response.  An  overall  smoothing 
of  the  beam  response  is  also  observed  for  this  problem. 
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Figure  5.  A  comparison  of  the  time-dependent  vertical  response  wo(L/2,t) 
(with  units  of  mm)  of  hinged-hinged  beams  due  to  a  suddenly  applied  trans¬ 
verse  load  q(t).  Results  are  for  both  viscoelastic  as  well  as  elastic  beams. 


time,  t 

Figure  6.  A  comparison  of  the  time-dependent  vertical  response  wo(L/2,t) 
(with  units  of  mm)  of  hinged-hinged  viscoelastic  and  elastic  beams  due  to  a 
periodic  concentrated  load  F(t)  (where  r)  =  2.744  x  108  N-sec./m2). 
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5.  Concluding  remarks 


In  this  paper  we  have  presented  an  efficient  finite  element  framework  for  the  numerical 
simulation  of  the  quasi-static  and  fully  transient  mechanical  response  of  viscoelastic  beam 
structures  based  on  the  kinematic  assumptions  of  the  third-order  Reddy  beam  theory.  The 
present  formulation  is  applicable  for  use  in  the  analysis  of  structures  undergoing  large  trans¬ 
verse  displacements,  moderate  rotations  and  small  strains.  The  viscoelastic  constitutive 
equations  in  the  semi-discrete  finite  element  equations  were  temporally  discretized  through 
the  use  of  the  trapezoidal  rule  and  a  two-point  recurrence  formula,  resulting  in  an  efficient 
temporal  integration  procedure.  The  finite  element  formulation  has  been  successfully  im¬ 
plemented  numerically  using  standard  low-order  reduced  integration  based  finite  element 
technology  as  well  as  a  family  of  elements  constructed  using  high  polynomial-order  Lagrange 
and  Hermite  interpolation  functions.  The  finite  element  procedures  have  been  successfully 
employed  in  the  numerical  simulation  of  the  mechanical  response  of  both  quasi-static  and 
fully  transient  viscoelastic  beams. 
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Appendix  A:  Semi-Discrete  Finite  Element  Equation  Components 

The  components  of  the  partitioned  element  coefficient  matrices  and  vectors  appearing  in  the 
semi-discretized  finite  element  formulation,  given  in  Eq.  (14),  may  be  determined  as 

rxb 

Mf1  =  /  Io^^dx  (Ala) 

J  Xa 

Mj2  =  M]l  =  0 
=  Mf l  =  0 
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dx 


X=Xb 
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In  the  above  equations  we  have  made  extensive  use  of  the  following  constants 


As  =  D0  —  2D2C2  +  D4C2 
L4  =  Cl(D4  —  DqCi) 

M2  —  D2  —  2D4Ci  +  DqC4 


(A4a) 

(A4b) 


(A4c) 


(A5) 


Appendix  B:  History  Variables  and  Fully-Discrete  Finite  Element 

Equation  Components 

The  history  variables  in  Eq.  (27)  may  be  expressed  explicitly  at  tk+ 1  as 
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Likewise,  the  components  of  the  fully  discretized  finite  element  coefficient  matrices  and 
vectors,  given  in  Eq.  (28),  may  be  evaluated  at  the  current  time  step  tk+ 1  as 
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The  quantities  (tfe+i)  are  expanded  as 
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where  n {  =  1,  n2  =  3  and  n3  =  2,  such  that  the  components  ]Qf  may  be  dehned  as 
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Appendix  C:  Element  Tangent  Stiffness  Matrix  Components 

The  components  of  the  element  tangent  stiffness  coefhcent  matrix,  given  in  Eq.  (31),  may 
be  determined  as 
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(Cla) 


Clearly  the  tangent  stiffness  coefficient  matrix  is  symmetric. 
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